function [interApprox, interApproxAV, ...
   eulerApprox, eulerApproxAV, ...
   interDHedgeApprox, interDHedgeApproxAV, ...
   eulerDHedgeApprox, eulerDHedgeApproxAV] = test2MonteCarlo(steps)

% Set the option parameters.
interestRate    = 0.06;
priceZero       = 100;
sigma           = 0.2;
T               = 1;
deltat          = 1 / 200;
N               = T / deltat;
type            = 'call';
strike          = 99;
epsilon         = 0.1;

iter = 0;
for M = steps
    % Keep track at which position in the result vector the current app-
    % roximation should be stored. This is necessary since the number of
    % steps do not need to be continuous.
   iter = iter + 1;
    % Initalize the Random Number Generator (RNG).
   randn('state', 1);
   for j = 1 : M
       % Generate a random walk phi_t ~ N(0, 1), t <= T
       phi     = randn(N);
       % Generate the random variable for the anti variate.
       phiav   = -phi;
       % Compute S_{1 <= t <= T}
       inter      = InitintegralMonteCarlo(interestRate, deltat, T, ...
            sigma, phi, priceZero);
       interAV    = InitintegralMonteCarlo(interestRate, deltat, T, ...
            sigma, phiav, priceZero);

      
       % We compute the delta hedge on the following way
       % Delta = (V(S + epsilon) - V(S)) / epsilon
       % this way we only need one extra path.
       interDH      = InitintegralMonteCarlo(interestRate, deltat, T, ...
            sigma, phi, priceZero + epsilon);
       interDHAV    = InitintegralMonteCarlo(interestRate, deltat, T, ...
            sigma, phiav, priceZero + epsilon);



       euler     = InitrecursiveMonteCarlo(interestRate, deltat, T, ...
            sigma, phi, priceZero);
       eulerAV   = InitrecursiveMonteCarlo(interestRate, deltat, T, ...
            sigma, phiav, priceZero);
       eulerDH   = InitrecursiveMonteCarlo(interestRate, deltat, T, ...
            sigma, phi, priceZero + epsilon);
       eulerDHAV = InitrecursiveMonteCarlo(interestRate, deltat, T, ...
            sigma, phiav, priceZero + epsilon);
       % Use the last value of the random walk S_{1 <= t <= T}, S_T to
       % compute the payoff of the option.
       if strcmp('call', type)
           % Set the values for the inter approximation.
           valueInt   = inter(N + 1) - strike;
           valueIntAV = interAV(N + 1) - strike;
           valueIntDH = interDH(N + 1) - strike;
           valueIntDHAV = interDHAV(N + 1) - strike;

           % Set the values for the euler approximation.
           valueEul   = euler(N + 1) - strike;
           valueEulAV = eulerAV(N + 1) - strike;
           valueEulDH = eulerDH(N + 1) - strike;
           valueEulDHAV = eulerDHAV(N + 1) - strike;

           % This factor is needed to take into account that the
           % deltahedge boundaries are inverted for a put when considering a
           % call.
           dhfactor = 1;
       elseif strcmp('put', type)
           % Set the values for the inter approximation.
           valueInt   = strike - inter(N + 1);
           valueIntAV = strike - interAV(N + 1);
           valueIntDH = strike - interDH(N + 1);
           valueIntDHAV = strike - interDHAV(N + 1);

           % Set the values for the euler approximation.
           valueEul   = strike - euler(N + 1);
           valueEulAV = strike - eulerAV(N + 1);
           valueEulDH = strike - eulerDH(N + 1);
           valueEulDHAV = strike - eulerDHAV(N + 1);
           dhfactor = -1;
       end
       
       % The values are set to the proper option so now everything can be
       % computed.
       % The regular pay off.
       payOffIntegral(j)       = max(valueInt, 0);
       % The pay off using AV.
       payOffIntegralAV(j)     = 0.5 * (payOffIntegral(j) + ...
         max(valueIntAV, 0));
       % The lowerbound of the delta hedge.
       payOffDHApproxUInt(j)   = max(valueIntDH , 0);
       % The lowerbound of the delta hedge with AV.
       payOffDHApproxUIntAV(j) = 0.5 * (payOffDHApproxUInt(j) + ...
         max(valueIntDHAV, 0));
      
       % The same computations as with the integral except now for the Euler.
       payOffEuler(j)          = max(valueEul, 0);
       payOffEulerAV(j)        = 0.5 * (payOffEuler(j) + max(valueEulAV, 0));
       payOffDHApproxUEul(j)   = max(valueEulDH , 0);
       payOffDHApproxUEulAV(j) = 0.5 * (payOffDHApproxUEul(j) + ...
         max(valueEulDHAV, 0));
   end
       
   % Store the approximated pay off of the option in a 'result' vector.
   interApprox(iter)       = exp(-interestRate * T) * sum(payOffIntegral) / M;
   interApproxAV(iter)     = exp(-interestRate * T) * sum(payOffIntegralAV) / M;
   % Compute the delta hedge
   interDHedgeApprox(iter)    = exp(-interestRate * T) * ...
      sum(payOffDHApproxUInt - payOffIntegral) / (epsilon * M);
   interDHedgeApproxAV(iter)  = exp(-interestRate * T) * ...
      sum(payOffDHApproxUIntAV - payOffIntegralAV) / (epsilon * M);

   % Store the approximated pay off of the option in a 'result' vector.
   eulerApprox(iter)       = exp(-interestRate * T) * sum(payOffEuler) / M;
   eulerApproxAV(iter)     = exp(-interestRate * T) * sum(payOffEulerAV) / M;
   % Compute the delta hedge
   eulerDHedgeApprox(iter)    = exp(-interestRate * T) * ... 
      sum(payOffDHApproxUEul - payOffEuler) / (epsilon * M);
   eulerDHedgeApproxAV(iter)  =  exp(-interestRate * T) * ... 
      sum(payOffDHApproxUEulAV - payOffEulerAV) / (epsilon * M);
end

% vim: expandtab :ts=4:tw=80

